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Abstract In this paper, we describe two computational methods for calculating the 
^V cumulative distribution function and the upper quantiles of the maximal difference 

between a Brownian bridge and its concave majorant. The first method has two dif- 
ferent variants that are both based on a Monte Carlo approach, whereas the second 
uses the Gaver-Stehfest (GS) algorithm for numerical inversion of Laplace transform. 
If the former method is straightforward to implement, it is very much outperformed 
by the GS algorithm, which provides a very accurate approximation of the cumulative 
distribution as well as its upper quantiles. Our numerical work has a direct application 
in statistics: the maximal difference between a Brownian bridge and its concave ma- 
jorant arises in connection with a nonparametric test for monotonicity of a density or 
regression curve on [0, 1]. Our results can be used to construct very accurate rejection 
region for this test at a given asymptotic level. 
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1 Introduction 

Consider the regression model Y^ = /o(ii) + £i, where tj = i/n for i — 1, • • • ,n, and 
conditionally on the regressors ij's, ei,---,e n are i.i.d. ~ (0, <t ) with < on < oo. 
Suppose that we are interested in knowing whether the true regression curve /o is 
nondecreasing on some sub-interval of [0, 1]. [6] considered the nonparametric test based 
on the maximum difference between the cumulative sum diagram of the observations 
and its concave majorant, multiplied by y/n and divided by any consistent estimator 
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The intuition behind this test is as follows: Under the null hypothesis that / is 
decreasing on [0, 1], the function L fo(t)dt, x £ [0, 1], is concave, and hence the cu- 
mulative sum diagram of the data must be "very close" to its concave majorant as 
n — > oo. [5] showed that asymptotically, the Type I error of the test attains its max- 
imum when /o is constant on [0, 1] . Then, under this least favorable case, the test 
statistic converges weakly to the maximum difference between a standard Brownian 
motion on [0, 1] starting at and its concave majorant. If (B(t),0 < t < 1), B(0) = 
denote a standard Brownian motion on [0, 1] starting at and B its concave majorant 
on [0, 1], then Durot's test statistic converges weakly to 

M = sup (S(t) - B(t)) . (1) 

te[0,i] v ' 

A similar testing problem for densities was considered by [TT]. The test can be based on 
the maximum difference between the empirical distribution and its concave majorant, 
multiplied by \fn. When the true density is uniform on [0, 1] , this maximum difference 
converges weakly to the distribution of M as in the regression setting above. 

As proved in Proposition 4 (iii) in [fQ, one interesting property of the distribution 
of M is that we can replace B in JT} by a standard Brownian bridge; i.e, the dis- 
tribution of M is also that of the maximum difference between a standard Brownian 
bridge and its concave majorant. Furthermore, the random variable M can be given 
under a more useful form. Let M 3 denote the maximum of a Brownian excursion and 
(A^3,lj M3.2, ' ■ ') an infinite sequence of independent random variables distributed as 
M3. If (Ui, U2, ■ ■ •) is an infinite sequence of independent uniform random variables on 
[0, 1] and (Li, L2, • • •) the corresponding uniform stick-breaking process; i.e., 

Li := U u L 2 := (1 - U X )U 2 , L 3 := (1 - Ux)(l - U 2 )U 3 ,- ■ ■ 

then [4] proved in Theorem 1 that 

M = raaxs/nM 3J . (2) 

j 

Absolute continuity of M is an immediate corollary of p]L Two other corollaries will 
follow from the same equality in distribution giving formulae for F, the cdf of M. If 
F 3 is the cdf of M 3 , then 



F(x) = E 



n^ 



V x > (3) 



The expectation in the formula above is taken with respect to the stick-breaking pro- 
cess, and F3 is known to be given by 

00 
F 3 (x) = 1 - 2 ^ {4n 2 x 2 - 1) exp(-2n 2 x 2 ), V x > 0, 

n=l 

see e.g. [TD] and [5]. 

A second formula, which follows from Proposition 7 and Theorem 8 of 0j, gives F 
as a function of the inverse of a Laplace transform. Let K m denote the modified Bessel 
function of the second kind and order meN, and G the function defined by 
00 
G{t) = JJ exp (-4 (2\FhnK\ (2Vlnt) - K (2\/2ntj\ ) , t > 0. (4) 



Then, 



F(x) = L~ 



G(Vt) 



t 



\x 2 ) 



x>0 



(5) 



where L~ h(z) denotes the value of the inverse of Laplace transform of h at z. 

We describe in Section 2 and 3 the implementation of two variants based on a 
Monte Carlo approach and a Gaver-Stehfest algorithm for approximating the inverse of 
Laplace transform. If the Monte Carlo (MC) methods are easy to implement, they both 
require a very large number of simulations in order to obtain the same precision as the 
deterministic Gaver-Stehfest (GS) algorithm. For x > 0.33, GS is able to approximate 
very accurately the cumulative distribution of M at x using a multiple precision library. 
For values x below what it seems to be a cut-off point for both methods, it is difficult 
to get a precise approximation for the distribution of M. This problem does not affect 
the calculation of the upper quantiles which is one of the main motivations of the work. 
Although the MC approach is not as efficient as the GS algorithm, it seemed natural 
to describe it in the sequel. We only report the numerical results of the GS algorithm, 
however. Tables __| __| and __| below give approximated values of the distribution function 
of M on a grid of real numbers x such that 0.33 < x < 2.54 with a regular mesh 
equal to 0.01. The approximation was performed with a precision ensuring up to 60 
significant digits. A table of quantiles of order p £ {0.90, 0.91, ■ • ■ ,0.99} is given as 
well. This table can be compared to the Monte Carlo approximated quantiles obtained 
by [6]. All the code used in the numerical computations in this paper is available at 
http: //www. ceremade .dauphine.fr/~fadoua/bf2010_code/ 



2 Monte Carlo approach 

We consider two different MC-based algorithms. They have the advantage of being 
very easy to understand and implement. The first approach is straightforwardly based 
on the expression of the distribution function of M given in Q . Because of the infinite 
product in (J3]l , a first approximation due to the truncation of the product is introduced. 
Control of the error due to this approximation is important in order to obtain a good 
theoretical estimator. Let J > be some finite integer and consider the problem of 
estimating 

Fj(x) = E f[F 3 \ -4=| . ., .(j. 

_7'=1 

For x > and a given e G (0, 1/4), the following lemma gives a lower bound for J 
so that 



< Fj(x) - F{x) < 2e. 



(6) 



Lemma 1 The approximation error satisfies {_]) if 

-log(x 2 6 2 /2) 



J>Jo 



log(2) 



+ 1. 



(7) 



Proof. See Appendix. 

(c) (c) (c) 

For x > and a given e > 0, we draw C independent copies (L\ ,L^ , J3 



■^) 



for c = 1, • • • , C to estimate Fj(x) where J = Jo as given in Lemma [S] The resulting 
Monte Carlo estimator is 

1 C i 

^/,c(*)=£En> 

The computation of the distribution function F3 imposes yet another approximation 
due to the fact that it is denned through an infinite series. The number of terms in 
the approximating finite sum needs to be larger for smaller values of x. Now by the 
Central Limit Theorem, we have 




fC(F, lc (x)-Fj(x))^ d M(0,aj) 



with 



= Var 



Ip3 



Let z 1 _ a / 2 be the (1 — a/2)- quantile of a standard normal for some small a G (0, 1) 
Then, for C large enough the event 



Fj{x) G 



Fj,c(x) 7 =- l -,Ft C (x) + 



Vc 



Vc 



occurs with probability » 1 — a. Combining both the deterministic and Monte Carlo 
approximations and noting that a j g [0, 1], it follows that 



F(x) G 



Fj, c (x)-2e 



■v/2 



Vc 



Fj,c(x) + 



z l-a/2 



Vc 



occurs with at least probability « 1 — a. Hence, to ensure an error of order e, the sample 
size C should be chosen of order [1/e J . Therefore, very large sample sizes are needed 
to get accurate results. To give an order of magnitude, Table [TJ shows several values 
of Jo and C corresponding to desired precision targets. All the values are computed 
for 0.33 < x < 2.54, where 0.33 appears to be the numerical limit of what we can 
compute without violating the basic properties of a distribution function. This point 
will be brought up again in the next section. Note that the main purpose of Table [T] is 
to give an idea about how Jo and C behave as functions of the precision. For instance, 
a precision of order 10 is useless if the goal is to compute an approximation of the 
value distribution function of M at 0.33 since it is of order 10 as found with the 
GS algorithm. 

We use the above MC approach to estimate the distribution function of M for 
0.33 < x < 2.54 as well as the upper quantiles. The algorithm is implemented in 
C. This method turns out to be very slow for large sample sizes. Moderate sample 
sizes (of order 10 ) do not give the desired accuracy for small x. The estimates of the 
distribution function for large x (of order 0.80 and above) as well as the upper quantiles 
match with those obtained by GS algorithm (see next section). 



In the same vein, one can consider a second variant of MC. It is mainly based on 
the following result due to Kennedy 1976 (see Corollary on page 372): 

M 3 = d sup B hl {t) - inf B hl (t) 
*6[0,1] *6[0,l] 

where B r is a Brownian bridge on length 1. Now using the well-known Donsker approx- 
imation, the distribution of M3 can be approximated for large N by the distribution 
of the random variable 



V N = sup VN(G N (t)-t)- inf VN(G N (t)-t) 
te[o,i] *e[o,i] 

where G^v is the uniform empirical process based on N independent uniform random 
variables Ui, ■ • ■ , t//v in [0, 1] . Using the fact that G^r is a constant function between 
the order statistics Cm) < ••• < Ui n \, it can be easily shown that 

Vjy = ViV < max I — — Uia\ ) — min I — -; U/;\ 

\\<i<N \N W J l<i<N \ N (l) 

Now the formula in ((2| yields the weak approximation 

M J,N — max ^LjVn 

1 < i < J 

where V^ , Vp, , • • ■ are independent random variables distributed as Vjsr , and J is a 
positive integer that should be chosen large enough to have the truncation error under 
control as done above. The distribution function of M can be estimated empirically 
by generating C independent random variables M\ L , M\ L, • • • , M\ 2, with the same 
distribution as Mjjy. If this second variant of MC has the drawback of adding another 
error due to the stochastic approximation of F3 by that of Vjv, it gives the possibility 
to generate samples with a distribution close to that of M for J, iV and C large enough. 
We will not pursue here the calculation of the approximation error as a function of J, 
N and C, which have to be very large to achieve high precision. 

The plot in Figure [5] shows an estimation of F using the first MC method with 
Jo = 100 and C = 10, 000. If the values are not accurate for small x, the plot gives 
nevertheless a good idea about the true shape of F. This is confirmed by the ap- 
proximation results we obtain with the numerical inversion of Laplace transform. The 
trajectory of 1000 independent random variables with the same distribution of Mj pf 
for J = 100 an iV = 10, 000 is shown in Figure [3] The sample was extracted from 
a larger one of size 10,000 with an empirical mean and standard deviation equal to 
0.9970 and 0.2475 respectively. 

If the MC approach gives a first idea of the support and shape of the distribution 
of M, it is not satisfactory in terms of efficiency and precision. As we show in the next 
section, the GS algorithm is a much better choice in both respects. 



3 Gaver-Stehfest algorithm 

The Gaver-Stehfest (GS) algorithm is one of several algorithms of numerical inversion 
of Laplace transform. For an excellent description of these algorithms, see pQ. The 
GS algorithm is different from other inversion procedures in that it involves only real 



numbers, but it also requires a very high numerical precision as we explain below (also 
see [J, p. 415). If g is the Laplace transform of some function / denned on R, then GS 
approximation of / is given by 

IK 

where K is an integer in N* and 

Under the assumption that the inverse of Laplace transform / has all its singularity 
points in (— oo,0] and that is infinitely differentiate on (0, oo), an extensive computa- 
tion study carried out by 3T has shown that 



/*■(*)-/(«) 



/(*) 



irr°- 8A ', t > 0. 



If the function / is bounded by 1 say, then the approximation in (JH| for well-behaved 
functions (in the sense given above) coincides with the truth up to significant 0.8K 
digits. Hence, the bigger K is, the better is the approximation. However, for large 
values of K, the binomial coefficients in £^ become extremely large and require high 
numerical precision. Such a facility is typically provided by a Multiple Precision (MP) 
numerical library or is built-in in some programming languages. 

For a given integer K > 0, let Fr denote the GS approximation of F. From the 
formula of F in ([5j) and (J8j , it is easily seen that 

2^ , 

F K (x) = Yl ^G{^k\o Z {2)x), x>0 (9) 

fe=l 

where G is the same function defined by the infinite product in Q. 

For x > and a given e > we approximate G by the product of the first N terms, 
where iV is a positive integer depending on x and e. Define 

N 
GjvW = n exp (-4 (2V2tnKi(2V2nt) - K a (2V2nt))) , t > 
3=1 

the truncated version of G. This truncation induces an additional error which we need 
to control. In fact, in computing the Gaver-Stehfest approximation of the distribution 
function F, we actually replace F^ in (|9j) by 

2 ^ e , 

F N ,k(x) = J2 T G n(V^o S (2) x), x>0. (10) 

fe=l 

The following shows that the error due to replacing F^(x) by Fnk(x) does not exceed 
a given threshold e > provided than JV is large enough. 



Lemma 2 For e > 0, we have \Fpf k( x ) ~ Fk( x )\ < £ tf N > Nq where 

N Q = - l \ In | 1 - ) + {2K + 1) \n(K) + 3K + 2 

+1. 

(11) 

Proof. See Appendix. 

From Lemma [2] it follows that 

\F NiK (x) - F(x)\ <e+\F K (x)-F(x)\. 

The second term in the left side is known to be of order 10 , and hence the 

approximation is of the same order if e is chosen to be o(10 ), and of order e if 

the latter dominates and N is chosen to be larger or equal than Nq given in (lll|) . 



We implement the multiple precision calculation of Fk in C++ using two open- 
source libraries for arbitrary precision computation: the GNU Multiple Precision Arith- 
metic Library (see [9]) and the Multiple Precision Floating-point Reliable Library 
(MPFR); see [7|. GMP is an optimized library written in C with assembly code for 
common inner loops. MPFR is built on top of GMP and adds support for common 
floating-point operations such as exp(s). 

To approximate the Bessel functions in Q, we use Bessel routines from the ALGLIB 
librarjU based on piecewise rational and Chebyshev polynomial approximations. We use 
a precision of 4000 bits to represent multiple precision floating-point numbers. However, 
the provided AGLIB Bessel approximations only guaranty a maximal error of order 
10 . As a proof-of-concept, we have also implemented the same algorithm using a 
much slower but more accurate numerical library in PythorjJ. For small values of x such 
as 0.30, 0.31, and 0.32, and unlike with the C library, we obtain results consistent with 
the monotonicity and positivity of a cumulative distribution function. For K = 60, 
N = 3200, the Python code gives the following approximations 9.8605317729e— 14 for 
x = 0.32, 1.10482969e-12 for x = 0.32 and 9.67030359e-12 for x = 0.33. 

Computing F^(x),x £ [0.33,2.54] takes about 6 hours (90 seconds per function 
evaluation) on a 2GHz single-processor machine. The computation is dominated by 
the evaluation of G in @. The coefficients £&, k — 1, • • • , 2K need to be computed only 
once. Tables [21 El and 0] give the approximated values of F on a grid starting at 0.33 
and ending at 2.54 with a regular mesh chosen to be equal 0.01. 

Finally, computing the upper quantiles of order is crucial when using the Kolo- 
mogorov type monotonicity test based on the maximal distance between the empirical 
cumulative sum diagram (resp. the empirical distribution) in the regression estimation 
setting (resp. the density estimation setting), see [6] and [11]. The GS algorithm can 
be easily used to approximate the upper quantiles of order p £ {0.90, 0.91, • ■ • , 0.99}. 

Note that these quantiles are between 1.33 and 1.72 (see Table[3]). For each quantile, 
we used a binary search and stopped when the difference between the GS approximation 



1 Available at http://www.alglib.net/specialfunctions/bessel.php 

2 Available at http://code.google.eom/p/mathpy/ 



of F at the point and the targeted probability falls below a given threshold (10 in 
the results we report). The results are shown in Table [5] This table is to be compared 
with the one published by [B] who obtained the quantiles for the same probabilities 
using a Monte Carlo approach. 



In this paper, Monte Carlo and a numerical inversion of the Laplace transform 
were used to estimate the distribution function and upper quantiles of M, the maximal 
difference between a Brownian motion on [0, 1] (or a Brownian bridge of length 1) and 
its concave majorant. This random variable determines the asymptotic critical region 
of a nonparametric test for monotonicity of a density or regression curve. We find the 
numerical inversion of Laplace transform, based here on the Gaver-Stehfest algorithm, 
to be much more accurate and faster than the Monte Carlo method. Numerical inversion 
of Laplace transform was then very well adapted to this problem. However, it would not 
have been possible to use such an efficient method if a Laplace transform representation 
of the distribution of M was not available, see H]. 

Finally, we would like to draw the reader's attention to the earlier computational 
work of [8] on Chernoff 's distribution. The latter appears as the limit distribution of 
the Grenander estimator; that is the Maximum Likelihood estimator of a decreasing 
density on (0,oo). In their work, [8] have also used a mathematical characterization 
of Chernoff's distribution. This allowed for a very efficient and fast approximation 
procedure which also outperformed Monte Carlo estimation. 



Table 1 Order of the lower bound Jo and sample size C. 



Precision 


Jo 


C 


10~ 5 


38 


10 -io 


10- 8 


57 


10 -16 


10 -io 


71 


10 -20 


10 -20 


138 


10 -40 




n 1 1 1 r 

0.0 0.2 0.4 0.6 0.8 1.0 






0.0 0.2 0.4 0.6 0. 



0.0 0.2 0.4 0.6 0.8 1.0 



Fig. 1 Four realizations of Brownian bridge and its concave majorant. The length of the dotted 
vertical segment equals M, the realization of the maximum difference between the Brownian 
bridge and its concave majorant. The Haar approximation was used to generate the Brownian 
bridge on a discrete partition of [0, 1] with a mesh equal to 2~ 12 . 



Appendix 

The following facts will be used in the proof of Lemma [T] 

Lemma A.l We have 

(i) For allj e N*, E(Lj) = 1/2-?. 
(U) For x > \/2, F 3 (x) > cxp(-l/x 2 ). 



Proof. The first identity can be proved recursively. For j — 1, we have E(L\) 
E{U\) = 1/2. Suppose that E(Li) = 1/2 1 for all i < j. It is easy to check that 



L j+1 = (1 - Ui)(l - U 2 ) ■ ■ ■ (1 - UrfUj+i = (1 - J2 L *) u 3+i 



10 



0.0 0.5 1.0 1.5 2.0 2.5 



Fig. 2 Plot of a Monte Carlo approximation of F based on a sample of size 10,000, with 
Jo = 100. 




Fig. 3 Plot of the trajectory of random sample of size 1000 of independent realizations of 
M JtN with J = 100 and N = 10,000. 
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Table 2 Approximated values of F obtained by the Gaver-Stehfest algorithm with K = 100, 
N = 3200. 



X 


F K 


X 


F K 


X 


F K 


0.33 


9.24257424322c-12 


0.61 


1.90415418636c-2 


0.89 


3.75627377529c-l 


0.34 


3.74465832649e-10 


0.62 


2.34012345624e-2 


0.90 


3.93227821925e-l 


0.35 


1.96857226601e-9 


0.63 


2.84117820291e-2 


0.91 


4.10762393335c-l 


0.36 


8.39648450672e-9 


0.64 


3.41073192977c-2 


0.92 


4.28198510241c-l 


0.37 


3.13734889037e-8 


0.65 


4.05156243499e-2 


0.93 


4.45505815829e-l 


0.38 


1.04444192578c-7 


0.66 


4.76576519986e-2 


0.94 


4.62656184022e-l 


0.39 


3.13540455103c-7 


0.67 


5.55472601039e-2 


0.95 


4.79623704507c-l 


0.40 


8.57933368778c-7 


0.68 


6.41911192791e-2 


0.96 


4.96384649825e-l 


0.41 


2.16022136991c-6 


0.69 


7.35887911016e-2 


0.97 


5.12917427370c-l 


0.42 


5.04741838646e-6 


0.70 


8.37329554403e-2 


0.98 


5.29202518946e-l 


0.43 


1.10248951989e-5 


0.71 


9.46097647495e-2 


0.99 


5.45222410267c-l 


0.44 


2.26592425775c-5 


0.72 


1.06199301858e-l 


1.00 


5.60961512572c-l 


0.45 


4.40745827182e-5 


0.73 


1.18476117686e-l 


1.01 


5.76406078288e-l 


0.46 


8.15505378035e-5 


0.74 


1.31409826206e-l 


1.02 


5.91544112452e-l 


0.47 


1.44191509163c-4 


0.75 


1.44965735553c-l 


1.03 


6.06365281379c-l 


0.48 


2.44619527193e-4 


0.76 


1.59105496305c-l 


1.04 


6.20860819881c-l 


0.49 


3.99630103610c-4 


0.77 


1.73787750362e-l 


1.05 


6.35023438140e-l 


0.50 


6.30744831947c-4 


0.78 


1.88968766384c-l 


1.06 


6.48847229167e-l 


0.51 


9.64597141746e-4 


0.79 


2.04603050312c-l 


1.07 


6.62327577632c-l 


0.52 


1.43309828337c-3 


0.80 


2.20643921928e-l 


1.08 


6.75461070708e-l 


0.53 


2.07334764606e-3 


0.81 


2.37044050697c-l 


1.09 


6.88245411417e-l 


0.54 


2.92727237058e-3 


0.82 


2.53755946188e-l 


1.10 


7.00679334912c-l 


0.55 


4.04100308041c-3 


0.83 


2.70732400181c-l 


1.11 


7.12762527962e-l 


0.56 


5.46401279238e-3 


0.84 


2.87926879124c-l 


1.12 


7.24495551883e-l 


0.57 


7.24806261908c-3 


0.85 


3.05293866924e-l 


1.13 


7.35879769041e-l 


0.58 


9.44600944518c-3 


0.86 


3.22789159067c-l 


1.14 


7.46917273011c-l 


0.59 


1.21105368429e-2 


0.87 


3.40370109977c-l 


1.15 


7.57610822418c-l 


0.60 


1.52928712783c-2 


0.88 


3.57995836076c-l 


1.16 


7.67963778457c-l 



By independence of (L\, ■ ■ ■ , Lj) and J/j+i, we can write 



3) aL ^ W J- 
i 



E(L J+l ) = E(1-J2 Li)/2 = (1 - Y, VW = 1/2 J+1 



i=l 



and the identity is proved for all j € N*. 

For the second inequality, we will use the fact that for a given a > 1 

(4a 2 t - 1) exp(-2a 2 t) < exp(-ai), for all t > 0. 

Consider the function 



(12) 



h(t) ~ (4a . t - 1) exp(-a(2a - l)t), t > 0. 

The study of variations of h shows that h is increasing on [0, (6a — l)/(2a — 1)] and 
decreasing on [(6a — l)/(2a — 1), oo) with with h(0) = — 1, linit^oo h(t) = and 



6a- 1 



2a - 1 / 2a - 1 



4a / 6a - 1 

— eXP (-^a- 



Now, the function 



log ft 



6a -1 
2a -1 



log 



4a 



2a -1 



6a -1 
4a 
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Table 3 Approximated values of F obtained by the Gaver-Stehfest algorithm. 



X 


F K 


X 


F K 


X 


F K 


1.17 


7. 7798004601 lc-1 


1.51 


9.62020662224e-l 


1.85 


9.96016083423e-l 


1.18 


7.87664018322e-l 


1.52 


9.64212812296e-l 


1.86 


9.96298535959e-l 


1.19 


7.97020525083e-l 


1.53 


9.66292598240e-l 


1.87 


9.96562371771e-l 


1.20 


8.06054783852e-l 


1.54 


9.68264838078e-l 


1.88 


9.96808708806e-l 


1.21 


8.14772354647e-l 


1.55 


9.70134203687e-l 


1.89 


9.97038605986e-l 


1.22 


8.23179097587e-l 


1.56 


9.71905221291e-l 


1.90 


9.97253065744e-l 


1.23 


8.31281133430e-l 


1.57 


9.73582272276c-l 


1.91 


9.97453036494e-l 


1.24 


8.39084806872c-l 


1.58 


9.75169594296e-l 


1.92 


9.97639415022e-l 


1.25 


8.46596652448e-l 


1.59 


9.76671282637e-l 


1.93 


9.97813048825e-l 


1.26 


8.53823362907e-l 


1.60 


9.78091291833e-l 


1.94 


9.97974738360e-l 


1.27 


8.60771759907e-l 


1.61 


9.79433437494e-l 


1.95 


9.98125239231c-l 


1.28 


8.67448766898e-l 


1.62 


9.80701398339e-l 


1.96 


9.98265264308e-l 


1.29 


8.73861384073e-l 


1.63 


9.81898718407e-l 


1.97 


9.98395485767e-l 


1.30 


8.80016665251e-l 


1.64 


9.83028809430e-l 


1.98 


9.98516537064e-l 


1.31 


8.85921696569e-l 


1.65 


9.84094953345e-l 


1.99 


9.98629014836e-l 


1.32 


8.91583576893e-l 


1.66 


9.85100304937e-l 


2.00 


9.98733480735e-l 


1.33 


8.97009399815c-l 


1.67 


9.86047894590e-l 


2.01 


9.98830463190e-l 


1.34 


9.02206237159e-l 


1.68 


9.86940631128e-l 


2.02 


9.98920459107e-l 


1.35 


9.07181123890e-l 


1.69 


9.87781304739e-l 


2.03 


9.99003935491e-l 


1.36 


9.11941044337e-l 


1.70 


9.88572589969e-l 


2.04 


9.99081331015e-l 


1.37 


9.16492919660e-l 


1.71 


9.89317048762e-l 


2.05 


9.99153057516e-l 


1.38 


9.20843596472e-l 


1.72 


9.90017133547e-l 


2.06 


9.99219501431e-l 


1.39 


9.24999836546e-l 


1.73 


9.90675190351e-l 


2.07 


9.99281025174e-l 


1.40 


9.28968307546e-l 


1.74 


9.91293461936e-l 


2.08 


9.99337968446e-l 


1.41 


9.32755574715e-l 


1.75 


9.91874090944e-l 


2.09 


9.99390649494e-l 


1.42 


9.36368093452e-l 


1.76 


9.92419123041e-l 


2.10 


9.99439366308(3-1 


1.43 


9.39812202742e-l 


1.77 


9.92930510053e-l 


2.11 


9.99484397768e-l 


1.44 


9.43094119365e-l 


1.78 


9.93410113095e-l 


2.12 


9.99526004728(3-1 


1.45 


9.46219932852e-l 


1.79 


9.93859705663e-l 


2.13 


9.99564431063(3-1 


1.46 


9.49195601129e-l 


1.80 


9.94280976712c-l 


2.14 


9.99599904647e-l 


1.47 


9.52026946811c-l 


1.81 


9.94675533688e-l 


2.15 


9.99632638303e-l 


1.48 


9.54719654107e-l 


1.82 


9.95044905522e-l 


2.16 


9.99662830687e-l 


1.49 


9.57279266289e-l 


1.83 


9.95390545586e-l 


2.17 


9.99690667143e-l 


1.50 


9.59711183695e-l 


1.84 


9.95713834588e-l 


2.18 


9.99716320502e-l 



is decreasing on [1, oo) with logft(l) = log(4) — 5/4 < 0, and hence h((6a— l)/(2a— 1)) < 
1. It follows that h(t) < 1 and the inequality in (|12[) is proved. 
It follows that 

oo oo 

1 - F 3 (x) = 2 ^2(4k 2 x 2 - 1) exp(-2k 2 x 2 ) < 2 V exp(-fcx 2 ) 



fc=l 



fc=l 



2exp(— x 



1 — exp(— x 2 ) 
To show that ^3(3;) > exp(— 1/x ) for all x > 2, it is enough to show that 

1 t^t- > exp(— 1/x ), for all x > \/2 

1 — exp(— x z ) 

or equivalently 

2cxp(-i) < (1 - exp(-t))(l - exp(-l/i)), for all t > 2. 
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Table 4 Approximated values of F obtained by the Gaver-Stehfcst algorithm. 



2.19 
2.20 
2.21 
2.22 
2.23 
2.24 
2.25 
2.26 
2.27 
2.28 
2.29 
2.30 
2.31 
2.32 
2.33 
2.34 
2.35 
2.36 
2.37 
2.38 
2.39 
2.40 
2.41 
2.42 
2.43 
2.44 
2.45 
2.46 
2.47 
2.48 
2.49 
2.50 
2.51 
2.52 
2.53 
2.54 



Fk{x) 



9.99739951848e-l 
9.99761711238e-l 
9.99781738392e-l 
9.99800163331c-l 
9.99817106996e-l 
9.99832681825e-l 
9.99846992298e-l 
9.99860135450e-l 
9.99872201360e-l 
9.99883273608e-l 
9.99893429703e-l 
9.99902741490e-l 
9.99911275534e-l 
9.99919093474c-l 
9.99926252361c-l 
9.99932804973c-l 
9.99938800114c-l 
9.99944282889e-l 
9.99949294965e-l 
9.99953874813c-l 
9.99958057939e-l 
9.99961877094e-l 
9.99965362474e-l 
9.99968541907e-l 
9.99971441026e-l 
9.99974083431c-l 
9.99976490838e-l 
9.99978683223e-l 
9.99980678951o-l 
9.99982494897e-l 
9.99984146562e-l 
9.99985648176c-l 
9.99987012798e-l 
9.99988252405e-l 
9.99989377977c-l 
9.99990399575e-l 



Table 5 Approximated upper quantilcs qi- a of order 1 — a. The approximation is based on 
the Gaver-Stehfcst algorithm with K = 100, N = 60. 



Q 

qi-a 



0.01 
1.71974853 



0.02 
1.61439819 



0.03 
1.54926391 



0.04 
1.50122253 



0.05 
1.46279052 



0.06 
1.43055908 



0.07 
1.40267791 



0.08 
1.37802490 



0.09 
1.35586822 



0.10 
1.33570159 



The preceding inequality can be proved as follows. Define the function 

k(t) ■- (exp(t) - 1)(1 - cxp(-l/i)), t > 2. 
We will show now that k(t) > 2 for all t > 2. For t > 2, we have 

1 \ / , ,.xl . exp(-l/t) 



k'(t) = exp(t) {l - (l + -) exp(-l/t)} + 



t 2 
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> ex 



p(t){l-(l + l)exp(-l/t)} 

exp(t)0(l/t) 



where 



<t)(z) = l-(l + z 2 )exp(-z), z£ [0,1/2]. 

It is easy to show that cj> is increasing on [0, 1/2] and hence 4>{ z ) > 0(0) = 1. It follows 
that the function k is increasing on [2, oo). Since fe(2) w 2.514 > 0, the inequality 
F 3 (x) > exp(-l/x 2 ), x > y/2 follows. D 

Proof of Lemma\T\ Define 



a/-i- n *(-s 

3=J+1 



We have 



0<^j(s)-F(a;) = E 



n* 



< E[Aj] 

= E[Ajl Aj < e ] +E[Ajl Aj>e ] 

<e + P(A,;>e). 



Let A j be the event 



Aj = {-kj < ^V 4 - fo r all j > J + l| 



and A j its complement. 
We can write 



P(Aj>e) = p( H F 3 



<l-e 



< 1 - e } n Aj )+ P(A 



OO 

/'i < n ^ 

jW+i 



< P( Yl exp^Lj/x 2 ) < 1 - e) + P(Aj), using Lemma A.l (i 

= p( f; L J >x 2 log(l/(l- e ))]+P(A c J ) 

/ OO v OO 

<P( £ L i>a ; 2 log(l/(l-e))) + £ P(L J>a; 2 /4). 
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Using Lemma A.l (ii) and the Chebyshev inequality, we get 

E.w+il/2 3 



p( J2 ^>^ 2 iog(i/(i- e ))) < 



z 2 log(l/(l- e )) 



2 J x 2 log(l/(l - e)) 



and 



j=J+i 



2 J x 2 



Hence, 



< Fj(x) -F(x) < e + 



+ 



2 J ir 2 log(l/(l-e)) 2 J 



To have this approximation error smaller than 2e, it suffices to take 



J> 



log(2) 

If e < 1/4 we can take 

and Lemma \fi\ is proved. 



lop 



J> 



+ log 



log(l/(l-e)) 



+ 4 



-log(xV/2) 
log(2) 



+ 1 



(13) 

□ 



Proof of Lemma The modified Bessel function of the second kind K n is known to 
converge to as x — > oo. Moreover we have 

and 



xK\{x) — Kq(x) < \[—^Jxexp(— x), \/x > 



see Lemma A. 2. For t > 0, define 



H(t) = 4 ^ (2V2ntK 1 (2V2nt) - ^(2^2^)] 



so that G(t) = exp(-H(t)). Also, for N G N* let 



ffiv(*) = 4 ^ ( 2V2ntK 1 (2V2nt) - K (2V2nt) 
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so that Gjv(i) = exp(— H]y(t)). We have 

< G N (t) - G(t) = exp(-H N (t)) - exp(-H(x)) < H(t) - H N {t) 

oo 
= 4 V^ (ncKi(nc) - K (nc)) 

- V 2 1 - exp(-c/2) V ' 

where c = 2y2t. 

Let us write again Fjv if f° r the Gaver-Stehfest approximation of the inverse of 
Laplace transform of Gjv- 

The corresponding approximation error due to truncating G is given by 

IK 

E NtK (x) = F K (x) - F NtK (x) = Y,T ( G (\/ feln ( 2 ) a; ) ~ G w ( v / fcln(2)x)) , x > 0. 



fc=i 



By (|14|) . we can write 



l^.^(x)h4./fgM exp(-a fc iV) 

' V 2 ri fc x ~ ex p(- a fe) 

A:— 1 



where a& = ^/21n(2)fe a;. 

Now, exp(-a fc ) < exp(-i/2ln(2)a:) and so (l-exp(-a fc ))~ 1 < (l-exp(-i/21n(2)a:)) _:L 
for fc = 1, • • • , 2A'. The coefficients fj. can be loosely bounded using the following upper 
bounds for binomial coefficients 

B ^<f«r, and M<^. 

For A: = f , • ■ • , 2K , we have 

kAK , T , N K 



i&i^i E "' +1 ^ (*w fc 



J=L(fc+l)/2J 

kAK 

k+l (T s s K(r, n \k 

K\ 

J=L(fc+l)/2J 



-TFT E fc^W^-i, 



so that 



< ^|^ +2 (2e 2 )^ 



^M<_L K 2K+1 (2e 2 ) K 
fc=l 



Hence, if we impose that \Ej^jq{x)\ < e, then it is enough to choose N such that 

N > - l \ In | l - ) + (2K + 1) ln{K) + 3K + 2 } . D 
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Lemma A. 2 For all x > 0, we have 

xK\(x) — Kq{x) < * / — v^exp(— x). 

Proof. Let us recall some well-known facts about modified Bessel functions of the 
second kind. 

K lM*) = J%^T-> forall^GC* (15) 



/k exp(— z) . . _. 
7= — -, as z -> oo and n£n 
2 y/z 



K n (z) sa a I \= — ', as \z\ — > oo and n G N (16) 

V 2 y/z 

lim K n (x) = oo for all 77 GN (17) 

a:\iO 

A"i(z) » i as |z| \0 (18) 

(/J?„(z))' = -z n iY„-i(z), for all z£ C and 77 £ Z (19) 

K' n (z) = --#„(*) - K n+1 (z) for all 2 G C* and 77 G Z (20) 

z 

K v (x) < K v ,(x) for all x > and 7/ < 1/ G K. (21) 

see e.g. Abramowitz and Stegun 1964. Note first that by (|15p . the inequality stated in 
the lemma is equivalent to 

xK\(x) — Kq(x) < xK 1 / 2 (x), x > 0. 

From JTSJl, HJiJ) and ([T7)l. it follows that 

lim (a;i4Ti(:r) — Ko(x) — x^ ; 2 (a:)) = —00 and 
x\0 

lim (as-K^a;) - K Q (x) - xK 1/2 (x)) = 0. 

x— >-oo ' 

Let us write VK 1 ) = a^i(aO — Kq(x) — xK 1 / 2 (x),x > 0. Suppose now that there exists 
x > such that ip(x) > 0. This would imply that there exists y > such that tp(y) > 
and i//(y) = 0. Now, using JTSJ and (0Oj it follows that 

iP'(x) = -iifo(i) + ^l(x) + (x - -) K 1/2 (x), x>0. 

Hence, y satisfies 

yKi(y) - K (y) > yK 1/2 (y) and 

Ki(y)=(^-y)K 1/2 (y)+yK (y). 

It follows that 

(y 2 -l)K a (y)>y(y+^jK 1/2 ( y ). 

Since Kq(x) > and K 1 / 2 (x) > for all x > 0, we must have j/ > 1. But if y > 1, 
then the previous inequality implies 

tfofo) > V( y_ 1 f ) jf i/2(y) > ^i/a(y) 

which is impossible by (|21[) . 
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